library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(readxl)
library(plotly)
## Warning: package 'plotly' was built under R version 3.5.1
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.5.1
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
Al tener una coleccion de datos observable, y estudiarlo a traves de correlaciones, podemos encontrar que puede tener 1 o mas parametros que rigen el comportamiento de los datos. La herramienta que nos permite el analisis de los parametros y como estos impactan el comportamiento de los datos es Hiden Markiv models o H.M.M.
Los H.M.M son modelos de distribucion probabilisticos de proposito general para una o más variables. Y es especifico para casos de sieres de tiempo, categoricas y continuas.
se debe evaluar cual es la distribucion que mejor se apega a los datos, y para el proposito de este documento (estudio de terremotos) se utilizara la distribución de poisson.
A continuacion un de la cuenta de sismos por anio, estos datos nos permitiran ejemplificar las mezclas y como la necesidad de los datos para metricos nos insta a utilizar H.M.M.
| Año | Sismos por año |
|---|---|
| 1900 | 13 |
| 1901 | 14 |
| 1902 | 10 |
| 1903 | 9 |
| 1904 | 13 |
| 1905 | 18 |
| 1906 | 24 |
| 1907 | 18 |
| 1908 | 11 |
| 1909 | 19 |
| 1910 | 28 |
| 1911 | 18 |
| 1912 | 14 |
| 1913 | 18 |
| 1914 | 17 |
| 1915 | 17 |
| 1916 | 25 |
| 1917 | 18 |
| 1918 | 17 |
| 1919 | 13 |
| 1920 | 8 |
| 1921 | 12 |
| 1922 | 14 |
| 1923 | 22 |
| 1924 | 17 |
| 1925 | 15 |
| 1926 | 16 |
| 1927 | 17 |
| 1928 | 16 |
| 1929 | 16 |
| 1930 | 8 |
| 1931 | 23 |
| 1932 | 12 |
| 1933 | 13 |
| 1934 | 20 |
| 1935 | 20 |
| 1936 | 16 |
| 1937 | 20 |
| 1938 | 22 |
| 1939 | 17 |
| 1940 | 18 |
| 1941 | 20 |
| 1942 | 17 |
| 1943 | 32 |
| 1944 | 23 |
| 1945 | 14 |
| 1946 | 23 |
| 1947 | 14 |
| 1948 | 17 |
| 1949 | 16 |
| 1950 | 24 |
| 1951 | 13 |
| 1952 | 12 |
| 1953 | 15 |
| 1954 | 9 |
| 1955 | 12 |
| 1956 | 10 |
| 1957 | 23 |
| 1958 | 9 |
| 1959 | 12 |
| 1960 | 14 |
| 1961 | 11 |
| 1962 | 9 |
| 1963 | 17 |
| 1964 | 12 |
| 1965 | 13 |
| 1966 | 11 |
| 1967 | 10 |
| 1968 | 23 |
| 1969 | 15 |
| 1970 | 21 |
| 1971 | 18 |
| 1972 | 15 |
| 1973 | 10 |
| 1974 | 16 |
| 1975 | 14 |
| 1976 | 17 |
| 1977 | 14 |
| 1978 | 17 |
| 1979 | 14 |
| 1980 | 13 |
| 1981 | 11 |
| 1982 | 10 |
| 1983 | 15 |
| 1984 | 13 |
| 1985 | 13 |
| 1986 | 11 |
| 1987 | 14 |
| 1988 | 10 |
| 1989 | 6 |
| 1990 | 19 |
| 1991 | 13 |
| 1992 | 13 |
| 1993 | 12 |
| 1994 | 17 |
| 1995 | 20 |
| 1996 | 14 |
| 1997 | 17 |
| 1998 | 12 |
| 1999 | 18 |
| 2000 | 12 |
| 2001 | 14 |
| 2002 | 13 |
| 2003 | 14 |
| 2004 | 15 |
| 2005 | 11 |
| 2006 | 11 |
| 2007 | 14 |
library(dplyr)
library(plotly)
datos <- data.frame(anio =c(1900,1901,1902,1903,1904,1905,1906,1907,1908,1909,1910,1911,1912,1913,1914,1915,1916,1917,1918,1919,1920,1921,1922,1923,1924,1925,1926,1927,1928,1929,1930,1931,1932,1933,1934,1935,1936,1937,1938,1939,1940,1941,1942,1943,1944,1945,1946,1947,1948,1949,1950,1951,1952,1953,1954,1955,1956,1957,1958,1959,1960,1961,1962,1963,1964,1965,1966,1967,1968,1969,1970,1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983,1984,
1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007), frecuencia = c(13,14,10,9,13,18,24,18,11,19,28,18,14,18,17,17,25,18,17,13,8,12,14,22,17,15,16,17,16,16,8,23,12,13,20,20,16,20,22,17,18,20,17,32,23,14,23,14,17,16,24,13,12,15,9,12,10,23,9,12,14,11,9,17,12,13,11,10,23,15,21,18,15,10,16,14,17,14,17,14,13,11,10,15,13,13,11,14,10,6,19,13,13,12,17,20,14,17,12,18,12,14,13,14,15,11,11,14))
plot_ly(data = datos, x=~anio, y=~frecuencia, type = 'scatter', mode = 'lines+markers')
plot_ly(data = datos, x=~anio, y=~frecuencia, type = 'bar')
Al utilizar la distribucion de poisson \(p(x) = \frac{exp(-\lambda)\lambda^x}{x!}\) tenemos que la media es de 15, pero segun poisson la media y la varianza son iguales, para este ejemplo la varianza es de 20.33, no coinciden con la distribucion de Pisson, sin embargo el fenomeno satisface las condiciones para utilizar la distribucion de Poisson. Se plantea la hipostesis, que el resultado esta siendo afectado por 1 de 2 distribuciones de Poisson con medias \(\lambda_1, \lambda_2\).
La seleccion de que media impacta el proceso esta dada por un mecanismo aleatorio, suponemos que la media \(\lambda_1\) es seleccionada con probabilidad \(\delta_1\) y \(\lambda_2\) con probabilidad \(\delta_2= 1-\delta_1\).
Vamos a definir un mecanismo C para definir las mezclas donde:
\(C=\left\{ \begin{array}{lcc} 1 & con & probabilidad &\delta_1 \\ \\ 2 & con & probabilidad &\delta_2 \\ \end{array} \right.\)
La funcion de probabilidad para \(X\) esta dada por \[\sum_{i=1}{\delta_i p_i}\]
Para la estimacion de los parametros un de una distribucion mixta esta dada por la maximo parentesco (Maximum Likeleyhood) la cual es la misma es para casos discretos o continuos. \[L(\theta_1,...,\theta_m,\delta_1,...,\delta_m)= \prod_{j=1}^{n}\sum_{i=1}^{m}\delta_i p_i(x_j,\theta_i)\]. Los \(\theta\) son los vectores de los parametros de la distribucion y los \(\delta\) son los parametros de mezcla, los cuales suman 1; los \(x_i,...,x_n\) son las observaciones.
Para el caso de m=2 (2 componentes de la distribucion de Poisson) con medias \(\lambda_1\) y \(\lambda_2\), y sean \(\delta_1\) y \(\delta_2=1-\delta_1\). Tenemos entonces que la distribucion p es: \[p(x)=\delta_1 \frac{\lambda_i^{x_i} exp(-\lambda_i)}{x!}+(1-\delta_1)\frac{\lambda_2^{x_i} exp(-\lambda_2)}{x_i!}\]
Asumiendo que \(\delta_2 = 1-\delta_1\), los parametros a buscar son \(\lambda_1, \lambda_2,\delta_1\) y el maximo parentesco es:\[L(\lambda_1, \lambda_2, \delta_1 | x_1,...,x_n) = \prod_{i=1}^n (\delta_1 \frac{\lambda_1^{x_i}e^{-\lambda_1}}{x_i!} + (1-\delta_1) \frac{\lambda_2^{x_i}e^{-\lambda_2}}{x_i!})\]
mllk <- function(wpar,x){
zzz <- w2n(wpar)
-sum(log(outer(x,zzz$lambda,dpois)%*%zzz$delta))
}
n2w <- function(lambda,delta){
log(c(lambda, delta[-1]/(1-sum(delta[-1]))))
}
w2n <- function(wpar){
m<- (length(wpar)+1)/2
lambda<- exp(wpar[1:m])
delta <- exp(c(0,wpar[(m+1):(2*m-1)]))
return (list(lambda=lambda,delta=delta/sum(delta)))
}
likelihood <- function(parms){
return(prod(outer(x,params$lambda,dpois) %*% params$delta))
}
Veamos con \(\delta_1\) y \(\delta_2\)
x <- datos$frecuencia
wpar <- n2w(c(10,20),c(0.5,0.5))
params<-w2n(nlm(mllk,wpar,x)$estimate)
## Warning in nlm(mllk, wpar, x): NA/Inf replaced by maximum positive value
L = likelihood(params)
params <- c(params, likelihood=L, "log L"=-log(L))
as.data.frame(params)
## lambda delta likelihood log.L
## 1 22.59287 0.09397933 2.096784e-135 310.1086
## 2 14.59052 0.90602067 2.096784e-135 310.1086
Veamos con \(\delta_1\), \(\delta_2\) y \(\delta_3\)
x <- datos$frecuencia
wpar <- n2w(c(10,20,30),c(1,1,1)/3)
params<-w2n(nlm(mllk,wpar,x)$estimate)
## Warning in nlm(mllk, wpar, x): NA/Inf replaced by maximum positive value
L = likelihood(params)
params <- c(params, likelihood=L, "log L"=-log(L))
as.data.frame(params)
## lambda delta likelihood log.L
## 1 14.59060 0.31459256 2.096784e-135 310.1086
## 2 14.59050 0.59143019 2.096784e-135 310.1086
## 3 22.59295 0.09397725 2.096784e-135 310.1086
Veamos con \(\delta_1\), \(\delta_2\), \(\delta_3\) y \(\delta_4\)
x <- datos$frecuencia
wpar <- n2w(c(10,20,35,40),c(1,1,1,1)/4)
params<-w2n(nlm(mllk,wpar,x)$estimate)
L = likelihood(params)
params <- c(params, likelihood=L, "log L"=-log(L))
as.data.frame(params)
## lambda delta likelihood log.L
## 1 14.59058 0.25197174 2.096784e-135 310.1086
## 2 14.59051 0.47784478 2.096784e-135 310.1086
## 3 14.59061 0.17620764 2.096784e-135 310.1086
## 4 22.59299 0.09397583 2.096784e-135 310.1086
Como podemos apreciar a partir del \(\delta_2\) no hay ganancia, y el algoritmo esta descomponiendo la \(\delta_1\) para llenar las variables faltantes, por ello -log(L) da el mismo resultado.
Una mejor apreciacion de mejor ajuste es a traves de la esperanza matematica y la varianza \(E(X)=\sum_i\delta_i\sigma_i\) \(Var(X)=E(X^2)-(E(X))^2\) siendo \(E(X^2)=\sum_i\delta_i(\lambda_i+\lambda_i^2)\)
Una secuencia de variales aleatorias discretas \({C_t:t\in \mathbb{N}}\) son consideradas como una cadena de Markov si para todos los \(t\in \mathbb{N}\) se satisface \[Pr(C_{t+1}|C_t...C_1)=Pr(C_{t+1}|C_t)\]
Es decir, el estado actual depende solo y unicamente del estado anterior, o en otras palabras; no tiene memoria. Por razones de espacio podemos escribir \(C^{(t)}\) como la historio \((C_1,C_2,...,C_t)\) o sea \(Pr(C_{t+1}|C^{t})=Pr(C_{t+1}|C_t)\)
Las cantidades asociadas a los estados de las cadenas de Markov son las probailidades condicionadas tambien llamdas probabilidades de transicion. \(Pr(C_{s+t}=j|C_s=i)\)
Si las probabilidades no dependen de \(s\) entonces la matriz es llamada “homogenea”, de lo contrario es “no homogenea”. Asumiremos que las matrices son homogeneas a menos que se indique lo contrario. \(\gamma_{ij}=Pr(C_{s+t}=j|C_s=i)\)
La matriz \(\Gamma(t)\) es definida como la matriz con \((i,j)\) elemento \(\gamma_{i,j}(t)\)
Una matriz finita de Markov satisface la ecuacion Chapman-Kolmogorov: \[\Gamma(t+u)=\Gamma(t)+\Gamma(u)\] La ecuancion de Chapman-kolmogorov implica para todo \(t\in\mathbb{N}\), \[\Gamma(t)=\Gamma(1)^t\] sea esto la matriz de probabilidad de transicion de t-pasos es la t-esima potencia de \(\Gamma(1)\), la matriz de primer paso. La matriz \(\Gamma(1)\), es la abreviacion de \(\Gamma\) es una matriz cuadrada cuyas filas suman 1. \[\begin{pmatrix} \gamma_{11} & \dots & \gamma_{1n} \\ \vdots & \ddots & \vdots \\ \gamma_{m1} & \dots & \gamma_{mn}\\ \end{pmatrix}\] Donde \(m\) denota el numero de estados de la cadena de Markov. La sentencia de las filas sumando 1 puede ser escrito como \(\Gamma1' = 1'\) esto es, el vector columna \(1'\) es un eigenvector de \(\Gamma\) y corresponden a los eigenvalores 1. Nos referimos a \(\Gamma\) como la matriz de probabilides de transicion (t.p.m, por sus siglas en ingles)
La probabilidad no condicionada \(Pc(C_t=j)\) de una cadena de Markov, en un estado dado, en un momento t es de interes y es denotado por el vector fila:
\(u(t)=(Pr(C_t=1),...,Pr(C_t=m)), \in \mathbb{N}\) Es decir nos permite saber el estado actual de la matriz. Siendo \(u(1)\) la distribucion incial de la matriz, para deducir la distribucion en el tiempo \(t+1\) desde \(t\) solo multiplicamos por la matriz de probabilidad de transicion: \[u(t+1)=u(1)\Gamma\]
Una cadena de Markov con matriz de probabilidad de transicion \(\Gamma\) se dice que tiene una distribucion \(\delta\) (un vector fila con elementos no negativos) si \(\delta\Gamma=\delta\) y \(\delta1'=1\).
El primer requerimiento expresa la estacionalidad y el segundo que \(\delta\) es de echo una distribucion de probabilidad.
La matriz \[\begin{pmatrix} 1/3 & 1/3 & 1/3 \\ 2/2 & 0 & 1/3 \\ 1/2 & 1/2 & 0 \\ \end{pmatrix}\]
Tiene una distribucion estacionaria \(\delta=\frac{1}{32}(15,9,8)\). Siendo que \(u(t+1)=u(t)\Gamma\) una cadena de Markov empezando desde esta distribucion tendra la misma distribucion en cualquier punto en el tiempo, esto es una matriz de Markov estacionaria.
Una cadena de Markov irreducible (homogenea, discreta en tiempo, espacio de estados finitos) tiene una unica, estrictamente positiva distribucion estacionaria. El vector \(\delta\) con elementos no negativos es una distribucion estacionaria de la cadena de Markov con t.p.m \(\Gamma\), si y solo si: \[\gamma(I_m-\Gamma +U)=1\] Donde 1 es un vector de unos, \(I_m\) es la matriz identidad \(m\times m\) y \(U\) es la matriz \(m\times m\) de unos.
La matriz estacionaria puede ser encontrada removiendo una de las ecuaciones del sistema \(\delta\Gamma=\delta\) y reemplazandola por \(\sum_i\delta_i=1\)
Asumimos que los estados \(1,2,..,m\) son cuantitativos y no categoricos, la funcion de autocorrelacion de \({C_t}\), se asume que es estacionaria e irreducible, se puede obtener que:
Primero, v=(1,2,…,\(m\)) y V=diag(1,2,…,\(m\)) tenemos los enteros no negativos \(k\), \[Cov(C_t,C_t+k)0\delta V\Gamma^kv'-(\delta v')^2\]
Si \(\Gamma\) es diagonalizable, y los eigenvalores (diferentes a 1) se denotan por \(\omega_2,\omega_3,...,\omega_m\), entonces \(\Gamma\) puede ser escrito como: \[\Gamma=U\Omega U^-1\] Donde \(\Omega\) es la diag(\(1,\omega_2,\omega_23,...,\omega_m\)) y U corresponde a los eigenvectores de \(\Gamma\). Tenemos para los valores no negativos \(k\) \[Cov(C_t,C_{t+k})=\delta VU\Omega ^kU^{-1}v'-(\delta v')^2\] \[= a\Omega ^kb'-a_1b_1\] \[\sum_{i=2}^{m}a_ib_i\omega_i^k\] donde \(a=\delta VU\) y \(b'=U^{-1}v'\) \[\rho(k)\equiv Corr(C_t,C_{t+k})=\frac{\sum_{i=2}^{m}a_ib_i\omega_i^k}{\sum_{i=2}^{m}a_ib_i}\]
Una vez observamos una cadenas de Markov y queremos estimar las probabilidades de transicion, una solucion - pero no la unica - is encontrar las transiciones y estimar las probabilidades como frecuencias relativas.
Suponiendo entonces, que queremos estimar \(m^2-m\) parametros \(\gamma_{ij} (i\neq j)\) de un estado \(m\) de la cadena de Markov \({C_t}\) El “likelihood” condicionado de la primera observacion es: \[L=\prod_{i=1}^{m}\prod_{j=1}^{m}\gamma_{ij}^{f_{ij}}\] Siendo el Log de este: \[l=\sum_{i=1}^{m}\sum_{j=1}^{m}f_{ij} log \gamma_{ij}=\sum_{i=1}^{m}l_i\]
y maximizamos \(l\), maximizando cada \(l_i\) por separado. Sustituyendo \(1-\sum_{k\neq i} \gamma_{ik}\), diferenciando \(l_i\) respecto a la diagonal de la probabilidad de transición \(\gamma_{ij}\) e igualando a cero tenemos que: \[0=\frac{-f_{if}}{1-\sum_{k\neq i}\gamma_{ik}}+\frac{f_{ij}}{\gamma_{ij}}=\frac{f_{ii}}{\gamma_ii}+\frac{f_{ij}}{\gamma_{ij}}\] A menos que un denominador sea cero, \(f_{ij}\gamma_{ii}=f_{ii}\gamma_{ij}\) y \(\sum_{j=1}^{m}f_{ij}=f_{ii}\). Esto implica que en el maximo del “likelihood”, \[\gamma_{ii}=\frac{f_{ii}}{\sum_{j=1}^{m}f_{ij}}\text{ y } \gamma_{ij}=\frac{f_{ij}\gamma_{ii}}{f_{ii}}=\frac{f_{ii}}{\sum_{i=1}^{}m}f_ij\] El estimador \(\gamma=\frac{f_{in}}{\sum_{k=1}^{m}f_{ik}(i,j=1,...,m)}\), es el la probabilidad de transicion empirica, el estimador \(\Gamma\) satisface el requerimiento de la sumas de las filas se igual a 1
Suponiendo la cadena de Markov \({C_t}\) toma valores 0 y 1, y deseasmos estimar las probabilidades de transicion \(\gamma_{ij}\) de la secuencia de observaciones en que \(f_{ij}\) transiciones del estado \(i\) al estado \(j\) (\(i\),\(j\)=0,1) y \(f_{11}>0\) pero \(f_{00}=0\). entonces en las observaciones un cero es seguido por un uno. \(c=f_{10}+(1-C_1)\) y \(d=f_{11}\). Las probabilidades de transicions estan dadas por \[\gamma_{01}=1 \text{ y } \gamma_{10}=\frac{-(1+d)+((1+d)^2 +4c(c+d-1))^{\frac{1}{4}}}{2(c+d-1)}\]
En el caso donde las observaciones con estados finitos parecen no satisfacer las propiedades de Markov, una alternativa es utiliza cadenas de Markov de orden superior, satisfaciendo: \[Pr(C_t|C_{t-1},C_{t-2,...})=Pr(C_t|C_{t-1},...,C_{t-1})\] Aunque el modelo no es una cadena de Markov per se, podemos redefinir el modelo de forma que el resultante lo sea. Si dejamos \(Y_t=(C_{t-l+1},C_{t-l+2},...,C_t)\) entonces \({Y_t}\) es una cadena de Markov de primer orden en \(M^t\) en el espacio de estados de \({C_t}\).
Una cadena de Markov de segundo orden, es estacionaria si se caracteriza por las probabilidades de transicion: \[\gamma(i,l,k)=Pr(C-t=k|C_{t-1},C_{t_2}=i)\] y tiene distribucion estacionaria bivariable \(u(j,k)=Pr(C_{t-1}=j,C_t=k)\) satisfaciendo \[u(j,k)=\sum_{i=1}^{m}u(i,j)\gamma(i,j,k) \text{ y } \sum_{j=1}^{m}\sum_{k=1}^{m}u(j,k)=1\] El uso general de las cadenas de Markov de orden superior incrementa el numero e parametros del modelo; una cadena general de Markov de orden \(l\) en \(m\) estados tiene \(m^l(m-1)\) probabilidades de transicion independientes. Param \(m=2\) hay modelos equivalentes, pero \(M>2\) pueden ser representados en una variedad de patrones dependientes y estructuras de autocorrelacion. En ambos casos un incremento de 1 orden en la cadena de Markov requiere solo un parametro adicional.
Los modelos de Raftery, quien denomina modelos “distribuciones de transicion mixtas” estan definidos como los procesos \({C_t}\) que toman valores \(M={1,2,...,m}\) y satisfacen \[Pr(C_t=j_0|C_{t-1}=j_1,...,C_{t-l}=jl) = \sum_{i=1}^{l}\lambda_i q(j_i,j_0)\] Donde \(\sum_{i=1}{l}\lambda=1\) y \(Q=(q(j,k))\) es una matriz \(m\times m\) con elementos no negativos que sus filas suman 1; cuyo lado derecho de la ecuacion esta acotada entre ceo y uno para todos los \(j_0,j_1,...,j_l\in M\)
Considerando los visto en la serie de terremotos en el capitulo anterior, se puede apreciar que las observaciones no tienen un limite de cantidad. Debido a esto, utilizar la distribucion de Poisson es la decision natural. En el capitulo uno se considero como puede afectar la interpretacion de la serie como una mezcla de distribuciones de Poisson, Cada una con una media \(\lambda_m\) donde m se refiere a la serie de datos y cada distribucion tiene un posibilidad de \(\delta_m\) donde \[\sum_{ i = 1}^{\infty} \delta_i = 1\]
Un modelo como una mezcla independiente se considera que no representara correctamente el modelo de los terremotos debido a que por definicion estas mezcla no tienen dependecia alguna con los valores anteriores sin embargo de acuerdo con los resultados de la grafica ACF una dependencia muy baja entre las observaciones, pero por fines ilutrativos se obtendra el valor de la correlaciones.
muestra <- read_xlsx("EarthQuakeData.xlsx", sheet="Sheet2", skip = 2)[-109,] %>%
rename("x" = !!names(.[2])) %>%
select(x)
acf(muestra$x)
La cadena escondida de Markov es una mezcla particular debido a su dependencia. Donde toda observacion \(X_t\) tiene adicionalmente una categoria \(C_t\) la cual se indica a cual de las distribuciones pertenece mientras que se puede calcular las probabilidades de la siguiente manera: \[Pr(C_t|C_{t-1}), t = 2,3..\] \[Pr(X_t|X_{t-1}, C_{}) t = 2,3.\]
Para este modelo se deben considerar dos parametros, primero un parametro no observado del proceso, \(C_t\) que satisface todas las propiedades de Markov, y un segundo parametro dependiente del estado, \(X_t\) en el cual la distribución de \(X\) depende unicamente de estado actual, \(C\).
Si la cadena de Markov tiene m estados, entonces se le conoce como HMM de m-estados. En este caso estamos considerando que apesar que ya se calcularos los valores de \(\delta_i\), existe adicionalmente una matriz de transicion como: \[ \Gamma = \begin{bmatrix} \gamma_{11} & \gamma_{12} \\ \gamma_{21} & \gamma_{22} \\ \end{bmatrix}\]
En contraste, en un caso de mezclas independientes, la distribucion de \(C_t\), el estado en el tiempo \(t\), si depende de \(C_{t-1}\). Por lo tanto, para considerar como se pueden representar las distribuciones de maneras discretas y continuas se debe puede utilizar como base la siguiente expresion para el estado \(i = 1,2,3..\) \[p_i(x)= Pr(X_t=x | C_t = i) \]
\(p_i\) es la probabilidad de \(X_t\) si la cadena de Markov esta en el estado \(i\) en el momento \(t\). En caso de variables continuas es necesarios considerar que la probabilidad como una funcion de densidad. A estas distribuciones \(p_i\) se les consideran como distribuciones dependientes del estado, teniendose \(m\) estados.
Ocacionalmente se necesita la distribucion marginal de \(X_t\) y distribuciones marginales mas elevadas, tal como \((X_t, X_{t+k})\). Estos resultados deben ser derivados de los casos de la matrices homogenias pero que no son necesariamente estacionarias.
Para valores discretos observados \(X_t\), definiendo \(u_i(t) = Pr(C_t = i)\) para \(t = 1,2 ...T\) se tiene
\[Pr(X_t = x) = \sum_{i = 1}^{m} Pr(C_t = i) Pr(X_t = x| C_t = i) = \sum_{i = 1}^{m} u_i(t)p_i(x) \]
La expreción se puede escribir de forma matricial de la siguiente manera: [Pr(X_t = x) = (u_1(t)…u_m(t)) \[\begin{bmatrix} p_1(x) & & 0 \\ & \ddots & \\ 0 & & p_m(x) \\ \end{bmatrix}\] \[\begin{bmatrix} 1 \\ \vdots \\ 1 \\ \end{bmatrix}\]]
Donde \(P(x)\) es definida como la diagonal de los \(p_i(x)\). Considerando que \(u_i(t) = u(1)\Gamma^{t-1}\), se obtiene: \[Pr(X_t = x) = u(1)\Gamma^{t-1}P(x) 1^{`}\]
Pero si se consideran el cambio de estado como una distribución estacionaria entonces se puede reescribir teniendo en cuenta que: \(\delta\Gamma^{t-1} = \delta\) para todo \(t \in \mathbb{N}\) se obtiene:
\[Pr(X_t = x) = \delta P(x)1^{`}\]
Primero debemos observar que: \[E(X_t) = \sum_{i=1}^{m}E(X_t|C_t=i)Pr(C_t=i)\]
Que al considerarla como un caso estacionario seria: \[E(X_t) = \sum_{i=1}^{m}\delta_i E(X_t|C_t=i)\]
En base a la ecuacion de arriba y considerado las propiedades de los momentos y de las cadenas se obtienen los siguientes momentos:
\[E(X_t) = \delta_1 \lambda_1 + \delta_2 \lambda_2\] \[Var(X_t) = E(X_t) \delta_1 \delta_2 ( \lambda_1 +\lambda_2)^2 > E(X_t)\] \[Cov(X_t, X_{t+k}) = \delta_1 \delta_2 ( \lambda_1 +\lambda_2)^2(1-\gamma_{12}-\gamma_{21})^k , For k \in \mathbb{N} \]
Cabe mencionar que la formula resultante de la covarianza tiene una forma tal que cuando las medias son iguales \(\lambda_1 = \lambda_2\) la covarianza se vuelve cero.
El objetivo de la seccion desarrollar la formula del likelihood \(L_T\) de \(T\) observaciones consecuentes asumiendo que son el resultado de una HMM de m estados, Esta formula existe pero para obtenerse el resultado es necesario calcular al rededor de \(m^T\) terminimos de \(2T\) factores, aparentando requerir \(O(Tm^T)\). Sin embargo hace tiempo esta comprobado que se puede calcular con \(O(Tm^2)\) operaciones. Esto es posible mediante la estimacion de los parametros por medio de optimizacion numerica.
Consdierando que el likelihood de una HMM en general. Supongamos una secuencias de observaciones \(x_1, x_2, .... x_T\) generador por este tipo de modelo. Se busca la probabilidad \(L_T\) de observar esta secuencia, como calculado bajo el modelo m-estado HMM que inicialmente tiene una distribución \(\delta\) y \(\Gamma\) de la cadena de Markov.
El likelihood de una Cadena de Markov Escondida (HMM), se puede representar por la siguiente ecuacion: \[L_t = P_r (X^{(T)} = x^{(T)}) = \delta P(x_1) \Gamma P(x_2) ... \Gamma P(x_T)1'\]
Siendo \(\delta\) la distribución inicial y \(P(x)\) la matriz diagonal de \(m*m\) con densidad de probabilidad \(p_i(x)\)
\(L_T\) se puede calcular como \(L_T = \alpha_T1'\) tomando \(\alpha_1 = \delta P(x_1)\) y de forma recursiva:
\(\alpha_t = \alpha_{t-1} \Gamma P(x_t),\) para \(t=2,3,...,T\)
Si la cadena de Markov se asume como estacionaria \((\delta=\delta \Gamma)\), en ese caso \(\alpha_0 = \delta\), por lo tanto:
\(\alpha_t = \alpha_{t-1} \Gamma P(x_t),\) para \(t=1,2,...,T\)
El likelihood es un producto de matrices, no de escalares, es por ello que no se puede simplemente calcular el log del likelihood como la suma de los logs de todos sus factores. Se sugiere un metodo de calcular el likelihood que supone \(log(p+q)\), en donde \(p>q\),
\[log(p+q) = log(p)+log(1 + q/p) = log(p)+log(1+exp(log(q)-log(p)))\]
El logaritmo de \(L_T\) se calcula escalando el vector de probabilidades \(\alpha_t\), es por ello que el vector se escala en cada \(t\) entonces sus elementos se suman a 1, manteniendo la sumatoria de los logs de los factores escalados. Por lo tanto el \(log(L_T)\) se puede definir como:
\[log(L_T) = \sum_{t=1}^T log(w_t/w_{t-1}) = \sum_{t=1}^T log(\phi_{t-1}\Gamma P(x_t)1')\]
Para calcular el log-likelihood de las observaciones \(x_1, x_2,...,x_t\) dentro del modelo Poisson-HMM con, al menos 2 estados, con matriz de probabilidades de transicion \(\Gamma\), vector de medias de estados dependientes \(\lambda\), y distribucion inicial \(\delta\), se utiliza el siguiente codigo implementado en R.
alpha <- params$delta * dpois(x[1], params$lambda)
lscale <- log(sum(alpha))
alpha <- alpha/sum(alpha)
for (i in 2:T) {
#alpha <- alpha %*% Gamma*as.numeric(dpois(x[i], params$lambda))
lscale <- lscale+log(sum(alpha))
alpha <- alpha/sum(alpha)
}
lscale
## [1] -2.386378
Primero se debe realizar una reparametrizacion para evitar las restricciones. Los elementos de \(\Gamma\), los de \(\lambda\) y el vector de medias de estados dependientes en una Poisson-HMM estan sujetos a restricciones de no negatividad y de otro tipo.
Las restricciones se pueden clasificar en 2 grupos. El primer grupo de restricciones son aquellas que dependen del estado de la distrubucion que se haya elegido, por ejemplo, la probabilidad de exito de una distribucion binomial se encuentra entre 0 y 1.
En el caso Poisson-HMM las restricciones mas importantes son: 1. La media \(\lambda_i\) del estado dependiente de la distribucion debe ser no negativa, para \(i=1,...,m\) 2. Las filas de la matriz de probabilidades de transicion \(\Gamma\) debe agregarse a 1, y todos los parametros \(\gamma_{ij}\) deben ser no negativos.
Para hacer la transformacion de los parametros, vamos a definir \(\eta = log \lambda_i\), para \(i=1,...,m\) en donde \(\eta \in \mathbb{R}\)
Luego de obtener la maximizacion del likelihood sin restricciones, se pueden obtener los parametros con restricciones transformando nuevamente: \(\lambda_i = exp (\eta_i)\).
La reparametrizacion de la matriz \(\Gamma\) se hace de la siguiente manera: \[\sum_{j=1}^m \gamma_{ij} = 1 \hspace{1cm} (i=1,...,m)\]
\(\Gamma\) tiene \(m^2\) entradas pero unicamente \(m(m-1)\) parametros libres, por lo que las se puede transformar a numeros reales sin restricciones \(\tau_{ij}, i \neq j\).
Para el caso \(m=3\) se define la matriz \[T = \begin{pmatrix} - & \tau_{12} & \tau_{13} \\ \tau_{21} & - & \tau_{23} \\ \tau_{31} & \tau_{32} & - \\ \end{pmatrix}\]
en donde \(\tau{ij} \in \mathbb{R}\). De tal forma que \(g : \mathbb{R} \rightarrow \mathbb{R}^+\) sea estrictamente una funcion creciente, \[g(x) = e^x \hspace{0.5cm} o \hspace{0.5cm} g(x) = \Bigg \{ \begin{array}{} e^x & x \leq0 \\x+1 & x \geq0 \end{array}\]
Entonces: \[ v_{ij} = \Bigg \{ \begin{array}{} g(\tau_{ij}) & para & i \neq j \\1 & para & i = j \end{array}\]
y:
\(\gamma_{ij} = \frac {v_{ij}} {\sum_{k=1}^m v_{ik}} \hspace{0.5cm} para \hspace{0.5cm} i,j = 1,2,...,m\) y \(\Gamma = (\gamma_{ij})\)
Por lo tanto los parametros \(\eta_i\) y \(\tau_{ij}\) los llamaremos parametros de trabajo, mientras que los paramentros \(\lambda_i\) y \(\gamma_{ij}\) seran los parametros naturales.
Utilizando estas transformaciones para \(\Gamma\) y \(\lambda\) el calculo de los parametros de la maximizacion del likelihood se pueden calcular asi: 1. Maximizar \(L_T\) con respecto a los parametros de trabajo \(T = \{\tau_{ij}\}\) y \(\eta = (\eta_1,...,\eta_m)\) 2. Transformar los estimadores de los parametros de trabajo en los estimadores de los parametros naturales: \(T \rightarrow \Gamma\), \(\eta \rightarrow \lambda\).
A continuacion se presenta un ejemplo de la una funcion en codigo de R haciendo referencia a la transformacion de los parametros naturales de Poisson en parametros de trabajo y viceversa.
pois.HMM.pn2pw <- function(m,lambda,gamma,delta=NULL,stationary=TRUE)
{
tlambda <- log(lambda)
foo <- log(gamma/diag(gamma))
tgamma <- as.vector(foo[!diag(m)])
if(stationary) {tdelta <-NULL} else {tdelta<-log(delta[-1]/delta[1])}
parvect <- c(tlambda,tgamma,tdelta)
return(parvect)
}
pois.HMM.pw2pn <- function(m,parvect,stationary=TRUE)
{
lambda <- exp(parvect[1:m])
gamma <- diag(m)
gamma[!gamma] <- exp(parvect[(m+1):(m*m)])
gamma <- gamma/apply(gamma,1,sum)
if(stationary) {delta<-solve(t(diag(m)-gamma+1),rep(1,m))} else
{foo<-c(1,exp(parvect[(m*m+1):(m*m+m-1)]))
delta<-foo/sum(foo)}
return(list(lambda=lambda,gamma=gamma,delta=delta))
}